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ABSTRACT 

Wc construct atmosphere models for strongly magnetized neutron stars with surface 
fields B - 10^2 _ ;l01^ G and effective temperatures T^b ^ 10*^ - 10^ K. The at- 
mospheres directly determine the characteristics of thermal emission from isolated 
neutron stars, including radio pulsars, soft gamma-ray repeaters, and anomalous X- 
ray pulsars. In our models, the atmosphere is composed of pure hydrogen or helium 
and is assumed to be fully ionized. The radiative opacities include free-free absorption 
and scattering by both electrons and ions computed for the two photon polariza- 
tion modes in the magnetized electron-ion plasma. Since the radiation emerges from 
deep layers in the atmosphere with p 10^ g/cm"^, plasma effects can significantly 
modify the photon opacities by changing the properties of the polarization modes. 
In the case where the magnetic field and the surface normal are parallel, we solve 
the full, angle-dependent, coupled radiative transfer equations for both polarization 
modes. We also construct atmosphere models for general field orientations based on 
the diffusion approximation of the transport equations and compare the results with 
models based on full radiative transport. In general, the emergent thermal radiation 
exhibits significant deviation from blackbody, with harder spectra at high energies. 
The spectra also show a broad feature {AE/Ebi ^ 1) around the ion cyclotron reso- 
nance Ebi = 0.63 {Z/ A) {B/ 10^'^ G) keV, where Z and A are the atomic charge and 
atomic mass of the ion, respectively; this feature is particularly pronounced when 
Esi SfcToff- Detection of the resonance feature would provide a direct measurement 
of the surface magnetic fields on magnetars. 

Key words: magnetic fields ~ radiative transfer - stars: atmospheres - stars: magnetic 
fields - stars: neutron - X-rays: stars 



1 INTRODUCTION 

Neutron stars (NSs) are created in the core collapse and subsequent supernova explosion of massive stars and begin their 
lives at high temperatures, T ^ 10^^ K. As they cool over the next ~ 10^ — 10^ years, they act as sources of soft X-rays 
with surface temperatures <; 10^ K. The cooling history of the NS depends on poorly constrained interior physics such as the 
nuclear equation of state, superfluidity, and magnetic field (see, e.g., Tsuruta 1998; Yakovlev et al. 2001, for a review) Thus 
observations of surface emission from isolated NSs can provide invaluable information on the physical properties and evolution 
of NSs. In recent years, the spectra of several radio pulsars (e.g., PSR B1055-52, B0656-I-14, and Geminga) have been observed 
to possess a thermal component that can be attributed to the surface emission at temperatures in the range (2-10) X 10^ K, 
while a few other pulsars show thermal emission from hot spots at temperatures of a few times 10® K on the NS surface (e.g., 
for observations in the X-rays, see Becker & Triimper 1997; Becker 2000; in the extreme ultraviolet, Edelstein et al. 1995; 
Korpela & Bowyer 1998; in the optical band, Pavlov et al. 1997; Caraveo et al. 2000). In addition, several radio-quiet, isolated 
NSs, presumably accreting from the interstellar medium (see Treves et al. 2000) , have been detected in the X-ray and optical 
bands (e.g., Caraveo et al. 1996; Waher et al. 1996; Waher & Matthews 1997; Caraveo 1998; Mignani et al. 1998; Paerels 
et al. 2001). Recent observations by the Chandra X-ray Observatory have also revealed a number of compact central sources 
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in supernova remnants with spectra consistent with thermal emission from isolated NSsQ. Most interestingly, in the last few 
years, thermal radiation has been detected from a potentially new class of NSs, "magnetars," which possess superstrong 
magnetic fields {B <: lO" G) (Duncan & Thompson 1992; Thompson & Duncan 1993, 1995, 1996; see Thompson 2000 for a 
review) . 

Observations of soft gamma-ray repeaters (SGRs) and anomalous X-ray pulsars (AXPs) strongly suggest that these 
are NSs endowed with superstrong magnetic fields (e.g., Vasisht & Gotthelf 1997; Kouveliotou et al. 1998; Kouveliotou et 
al. 1999; see, e.g., Hurley 2000 for SGR review; Mereghetti 1999 for AXP review). According to the magnetar model, the 
X-ray luminosity from AXPs and the quiescent X-ray emission from SGRs are powered by the decay of a superstrong [B ~ 
10^* — 10^^ G) magnetic field (Thompson & Duncan 1996) and/or by the residual thermal energy (Heyl & Hernquist 1997a). 
Alternatively, it has been proposed that the AXP emission may originate from accretion onto a NS with B ~ lO^'^ G, e.g., from 
a fossil disk left over from the supernova explosion (see Alpar 2000; Chatterjee et al. 2000, and references therein). However, 
detailed timing studies of AXPs and SGRs (e.g., Thompson et al. 2000; Kaspi et al. 2001) and deep optical observations 
(e.g., HuUeman et al. 2000; Kaplan et al. 2001), as well as the many similarities between AXPs and SGRs, clearly favor the 
magnetar model for both AXPs and SGRs. Nevertheless, it should be emphasized that, like radio pulsars, the magnetic field 
strengths {B ~ 10^* - lO^'^ G) of magnetars are indirectly inferred and are based on the measurement of the spin periods 
(6 — 12 s) and period derivatives and on the assumption that the spin-down is due to the loss of angular momentum that 
is carried away by magnetic dipole radiation or Alfven waves (Thompson & Blaes 1998; Thompson et al. 2000). Thermal 
radiation has already been detected in four of the five AXPs (IE 1048.1-5937, IE 1845-0258, 4U0142-H61, RXS J1708-40) 
and in SGR 1900-1-14; fits to the spectra with blackbody or with crude atmosphere models indicate that the thermal X-rays 
can be attributed to magnetar surface emission at temperatures of (3 — 7) x lO'' K (see Mereghetti 1999; Perna et al. 2001). 
Clearly, detailed observational and theoretical studies of thermal emission can potentially reveal much about the physical 
conditions and the true nature of magnetars. Of particular importance are possible spectral features, such as the ion cyclotron 
line [see eq. (1.1)] studied in this paper (see also Zane et al. 2001), which, if detected, can provide a direct measurement of the 
surface magnetic field. In fact, it has been suggested that the lack of proton and electron cyclotron features in XMM-Newton 
observations of the isolated neutron star RX J0720. 4-3125 may lead to constraints on its surface magnetic field (Paerels et 
al. 2001). 

Thermal radiation from a magnetized NS is mediated by the thin atmospheric layer (with scale height ~ 0.1 — 10 cm and 
density ~ 0.1 — 100 g/cm^) that covers the stellar surface. The physical properties of the atmosphere, such as the chemical 
composition, equation of state, and especially the radiative opacities, directly determine the characteristics of the thermal 
emission. While the surface composition of the NS is unknown, a great simplification arises due to the efficient gravitational 
separation of light and heavy elements (with timescales on the order of seconds; see Alcock & lUarionov 1980 for a discussion 
of gravitational separation in white dwarfs). A pure H atmosphere is expected even if a small amount of fallback/accretion 
occurs after NS formation^. A pure He atmosphere results if H is completely burnt up, and a pure Fe atmosphere may be 
possible if no fallback/accretion occurs. 

Steady progress has been made over the years in modeling NS atmospheres. The first models of zero-field NS atmospheres 
were constructed by Romani (1987). Later works used improved opacity and equation of state data from the OPAL project for 
pure hydrogen, helium, and iron compositions (Rajagopal & Romani 1996; Zavlin et al. 1996). These works showed that the 
radiation spectra from light element (H or He), low-field {B ^ 10* G) atmospheres deviate significantly from blackbody spectra. 
Magnetic NS atmosphere modeling was first attempted by Miller (1992), who unfortunately adopted an incorrect polarization- 
averaging procedure for the radiative transport and included only bound-free opacities. So far the most comprehensive studies 
of magnetic NS atmospheres have focused on hydrogen and moderate field strengths of B ~ 10^^ — 10^** G (Shibanov et 
al. 1992; Pavlov et al. 1994; see Pavlov et al. 1995 for a review and Zane et al. 2000 for full, angle-dependent atmosphere 
models with accretion). These models correctly take into account the transport of different photon modes (in the diffusion 
approximation) through a mostly ionized medium. The opacities adopted in the models include free-free transitions and 
electron scattering, while bound-free opacities are treated in a highly approximate manner and bound-bound transitions are 
completely ignored. Because the strong magnetic field significantly increases the binding energies of atoms, molecules, and 
other bound states (see Lai 2001 for a recent review), these bound states have appreciable abundances in the atmosphere 
at low temperatures {T ^ 10^ K) (see Lai & Salpeter 1997; Potekhin et al. 1999). Therefore the models of Pavlov et al. 
are expected to be valid only for relatively high temperatures (T <; a few x 10® K), where hydrogen is almost completely 
ionized. The presence of a magnetic field gives rise to anisotropic and polarized emission from the atmosphere (see Pavlov & 
Zavlin 2000). The resulting magnetic spectra are softer at X-ray energies than the low-field spectra but are still harder than 
the blackbody spectra. Models of magnetic iron atmospheres (with B ~ 10^^ G) were studied by Rajagopal et al. (1997). 



* see G. Pavlov 2000, bttp://online.itp. ucsb.edu/onlinc/ncustars_cOO/pavlov/ 

^ As far as the radiation spectrum is concerned, the photons decouple from the matter at Thomson depth t ^ 10^ (see Fig. P); 



the 



total mass above this depth is less than 10 



16 
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Because of the complexity in the atomic physics and radiative transport, these Fe models are necessarily crude. Despite a 
number of shortcomings, these H and Fe atmosphere models have played a valuable role in assessing the observed spectra of 
several radio pulsars and radio-quiet NSs (e.g., it is now well recognized that fitting the thermal spectra by a blackbody in 
the X-rays tends to overestimate the effective temperature) and some useful constraints on NS properties have been obtained 
(e.g., Meyer et al. 1994; Pavlov et al. 1996; Zavlin et al. 1998). 

This paper is the first in a series where we systematically investigate the atmosphere and spectra of strongly magnetized 
NSs. Here we study the H and He atmospheres with relatively high effective temperatures {TcS ^ a few x 10*' K) so that bound 
atoms or molecules may be neglected. We focus on the superstrong field regime {B <: 10^'' G) and compare with results from 
models possessing weaker magnetic fields. We construct self-consistent atmosphere models in radiative equilibrium and solve 
the full, angle-dependent, coupled transport equations for both polarization modes. For comparisons, we also study models 
based on the diffusion approximation of the transport equations. Of particular interest is the effect of the ion cyclotron 
resonance, which occurs at energy 

^IS^°-'KI)(ioSg)^^^' ^'-'^ 

where Z and A are the charge number and mass number of the ion, respectively. As we show in this paper, the resonance gives 
rise to a broad {AE/E ~ 1) spectral feature around Esi, and this feature is particularly pronounced when Ebz ^ SkTcs- 
Obviously, the ion cyclotron feature can potentially provide important diagnostics for the physical properties of magnetars. 

While in the final stages of writing our paper for publication, we became aware of two recent papers that also explore 
H atmospheres in the superstrong field regime. Ozel (2001) included the effect of vacuum polarization (see Section |^ but 
neglected the ion effect on both opacities and polarization modes, while Zane et al. (2001) studied the proton cyclotron feature 
and adopted an approximate description of the polarization modes in the electron-proton plasma with vacuum polarization 
(see also Section P). 

In Section ^, we review the basic physics ingredients of our atmosphere models, including the radiative transfer equations, 
the equation of state, and the opacities. The numerical methods for both full transport and diffusion models and tests of 
the solution are discussed in Section |^. In Section ^, we present our numerical atmosphere models and spectra for different 
compositions (H or He), effective temperatures, and magnetic field strengths and geometries. Finally, Section ^ contains a 
discussion of the uncertainties in our models and possible future works. 



2 ATMOSPHERE MODEL: BASIC EQUATIONS AND PHYSICS INPUTS 

We consider an isolated NS with a plane-parallel atmosphere. This is justified since the atmospheric scale height H ^ 10 cm 
is much less than the NS radius i? « 10 km. The atmosphere is composed of pure hydrogen or helium (see Section |l|) . The 
temperature T is taken to be ^ 10^ K, and the atmosphere is assumed to be fully ionized (see Section ^ for discussion on 
ionization equilibrium). A uniform magnetic field B permeates the atmosphere, and the angle between the direction of the 
magnetic field B and the surface normal z is Qb- 



2.1 Radiative Transfer Equation 

In a magnetized plasma, there are two normal modes of propagation for electromagnetic waves. These are the extraordinary 
mode (X-mode, j ~ 1), which is mostly polarized perpendicular to the k-B plane, and the ordinary mode (O-mode, j — 2), 
which is mostly polarized parallel to the k-B plane, where k is the unit vector along the wave propagation direction (e.g., 
Meszaros 1992; see Section ^.5| ). When the normal modes are approximately orthogonal and their relative phase shift over a 
mean free path is large (i.e., large Faraday depolarization), the radiative transfer of the four Stokes parameters reduces to 
that of the two modes (see Gnedin & Pavlov 1974). The radiative transfer equation (RTE) for the specific intensity II of 
mode j is 

2 
i=l 

where Bi, = Bv{T) is the blackbody intensity, nf'^iyC) (in cm^/g) is the absorption opacity (corrected for stimulated emission) 
for mode j propagating along k, di<^'^{k!i —* kj)/df2 is the differential opacity for scattering from mode i in direction k' to 

(k) = ftf =(k) 



mode j in direction k, K'°'(k) — i<f"^{\C) + Kj'^(k) is the total opacity, and 



Kf(k) = ^K-(k) (2.2) 

1=1 

K-(k) = Jd^'^^^^-^f^- (2-3) 
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Defining the Thomson depth rhy dr — —pn^'dz, where kq' = YgaT /rrip — 0.4 Yc cm^/g is the Thomson electron scattering 
opacity in the absence of a magnetic field and Yc = Z/A is the electron fraction, equation ( |2.l[ ) becomes 



- ^-^^ - ^ > . " I'^i^ h (2.4) 



where ^ = k • z = cos 6*. In general, at a given frequency v and depth r, the intensity Ii(k) depends on /i and the azimuthal 
angle (j) when B is not along the surface normal, i.e., when Ob 0. In solving the RTE, it is useful to introduce the two 
quantities 

ii{k) = i[/^(k) + /^(-k)l (2.5) 
flik) = i[/^(k)-/^(-k)l, (2.6) 

which are evaluated only for fJ, > and are related to the mean specific intensity and energy flux (see Mihalas 1978 for the 
nonmagnetic case). As a result, the RTE becomes 

fm = M^l^ (2.7) 



dfii^ _ 2 9V.(k) _ <-(k) i?. 2 V /^i,^ rf^°^(k'»^ki) 



where dr^, — — pKj°'(k)d2 is the total optical depth and implies that the domain of integration is restricted to k' ■ £ = 
/i' > 0. In deriving these equations, we have used the reflection symmetry of the opacities, i.e., Kj°'(k) = Kj°'(— k) and 
dt^^Ck'i ^ kj)/dn = dK^'^ik'i ^ -kj)/dQ. = dK'^'i-k'i ~kj)/dQ.. The boundary conditions for the RTE are Ii{k) = for 
^ < at r,^ = and Ii{k) By 12 at ^ oo. In terms of ij, these are 

= ii atr. =0 (2.9) 

OTy 

f^^ + ^i = Jf^l^ + B.) atr.^oo. (2.10) 

OTy 2 V dry J 

To obtain the constraint imposed by radiative equilibrium, we integrate the RTE ( |2.l[ ) over the solid angle, which yields 
the zeroth-order moment of the transfer equation, 

2 

V-Fi(k) =p J dk>,f^{k)^-p J dk^f\k)Ii{k) + pY^ JdkJ dk'^^^-^:^^^^f^ll{k'), (2.11) 
where the specific energy flux of mode j is 

Fi = J dkkli{k) = 2 J dkkflik). (2.12) 

In radiative equilibrium, the total flux, F = X^j^^ / diyFi, satisfies V ■ F = 0, so that 

2 

V-F = 2p^ f du f dkKf=(k) ^~ii{k) =0. (2.13) 



2 

Also, for an atmosphere with an effective temperature Toff, the total flux is constant and satisfies 

2 „ „ 2 



F-^^jdyj dk^ilKk) = '2^JduJ dkfifiik) = assT, 
where asB is the Stefan-Boltzmann constant. 



off) 



(2.14) 



2.2 Special Case: Magnetic Field Parallel to Surface Normal 

In the special case when B is parallel to the surface normal z, the intensity Ii{k) is independent of (p, so that Ii{k) — Il[p). 
In the atmosphere models based on the full RTE presented in this paper, we consider only this special case. It is clear that 



the RTE (^.8) can be written in the familiar form 



2 dii 



t^'j^=ii-Si, (2.15) 
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where 



SIM 



47r 



'(m) 2 



2 

E 



dfi' 



(2.16) 



is the source function. 

We note that, in using equation ( |2.16| ) to calculate the source function, the integral over the initial photon direction fj,' 
must be carried out at every frequency i/, depth r, and final photon direction /j, and for each global iteration (temperature 
correction) as described in Section ^ We can reduce the computation time by assu ming that the differential scattering cross- 
section is approximately independent of the initial photon direction. Equation ( 2.16 ) then yields the approximate source 
function 



Slip) 



(m) 2 



+ 



Eh 



'(^) in 



(2.17) 



where the specific energy density is 



= - / dkllili). 



(2.18) 



We expect that this approximation introduces only a small error to the emergent spectra, and a comparison of the final results 



using the source functions given by equations (2.16) and (2.17) indeed show the differences are small, though this has not 
been studied extensively (see also Shibanov & Zavli n 1995). The numerical results presented in Section ^ make use of the 
approximate source function given by equation ( 2.17 ) . 



2.3 Diffusion Approximation 



In this paper, we shall also present atmosphere models for general Os based on the diffusion approximation of the RTE. To 
obtain the diffusion RTE, we multiply both sides of the RTE ( |2.l[ ) by — k/pKj°'(k) and integrate over the solid angle to find 
the specific flux. 



P4°*(k) 



k(k- V)/^(k), 



(2.19) 



where we have used the symmetry property of the opacities, so that J dk [k/«:*°*(k)][d«:°'^(k'i kj)/dS7] — and 
J dkk[Kj''''(k)/Kj°'(k)] = 0. In the diffusion limit, we approximate the specific intensity by ~ {c/4n)[ul + (3/c)k ■ Fj]. 

Substituting into equation (2.19) and using the fact that VuJ — [dwl / dz)z and that the opacity only depends on ps = 
k- B 



cos Sfl, we find 



d 



f (^/^' cosQbB + if- sin esBj 



(2.20) 



where B± is the unit vector perpendicular to B and lies in the plane of B and z (recall that 0s is the angle between B and z) 

are the angle- averaged mean free paths parallel and perpendicular 



and B± sin Os = z — (z ■ B)B. In equation (2.2C ), ij and if 
to the magnetic field, respectively, and they are given by 



I, 



It = 



d^i 



2 

Ms 



p'^r(/^s) 



dpi 



p4°*(ms)' 



(2.21) 



(2.22) 



The specific flux perpendicular to the stellar surface is -F^^z « —clj{dui/dz), where the averaged mean free path along z is 



Z" cos Qs + li sin 0s. 



(2.23) 



Note that, as a consequence of equation (2.20), a magnetic field can induce horizontal radiative flux; this implies that a 
horizontal temperature gradient may develop in the atmosphere. To properly address this issue, one needs to study the three- 
dimensional structure of the atmosphere; this is beyond the scope of this paper (see Section Substituting the approximate 
expressions of the specific intensity and flux into the zeroth-order moment equation (2.11) and defining the Thomson depth 
dr = ~pkq^ dz and the zero-field Thomson mean free path Iq = l/(pK™), we obtain 



d_ 

dr 



Ij dul 



Kf 



+ 



(2.24) 
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where — [A-k / c)Bu{T) is the blackbody specific energy density and 



47r 



dk (k) 



This equation has been previously derived by Kaminker et al. (1982, 1983). 



The boundary conditions for equation (2.24) are £ ■ FJ = ciiJ/2 at r = and ^ ?i^/2 at r cxj, i.e., 



Ij dui 
lo dr 



duj 
dr 



+ ui 



T 

1 f dul 



2 V dr 



at T = 



+ uZ 



at r ^ oo. 



Finally, in the diffusion approximation, the radiative equilibrium condition [see eq. (2.13)] simplifies to 

2 



y-abs / p 
U 



(2.25) 
(2.26) 
(2.27) 

(2.28) 
(2.29) 

(2.30) 



2.4 Hydrostatic Equilibrium and Equation of State 

In modeling the atmosphere, we must determine the density p(r) and temperature T{t) along with the radiation intensity. 
From hydrostatic equilibrium, the pressure P at Thomson depth r is given by 
9 , 



P(r) 



(2.31) 



where g = (GM/R^){1 - 2GM / Rc^)~^^'^ = 2.4 x lO^"* cm/s^ is the gravitational acceleration at the surface of a M=1.4 M©, 
_R=10 km neutron star. We assume the atmosphere is fully ionized. The equation of state can be written as 



P = P, + Pon + Pco 



(2.32) 



where Pc is the pressure due to electrons, which may or may not be degenerate, Pion is the pressure due to ions, which 
still behave as a classical ideal gas under the conditions that exist in the upper layers of the NS, and Pcoui is the Coulomb 
correction to the pressure. In strong magnetic fields, electrons are quantized into Landau levels, and the electron pressure is 
given by (e.g., Lai 2001) 



B 



v=0 



(2.33) 



where Ac = h/rricC is the electron Compton wavelength, Pq = rri^c'^ / {eh) = 4.414 x 10^^ G, is the spin degeneracy of the 
Landau level [go = 1 and g^ = 2 for v > 1), fic is the electron chemical potential, and F„{y) is the Fermi-Dirac integral 



Fniy) 







-y + 1 



-dx. 



When the conditions 



T Tb 

p < pB 



k 



1.34 X 10" Pi2 K 



7.05 X lO'' Y-^B%^ g7cm^ 



(2.34) 

(2.35) 
(2.36) 



where Pq = [hc/eBf^'^ and P12 = P/(10^^ G), are satisfied, as in the case of the atmospheres considered in this paper, only 
the ground Landau level (v = 0) is occupied. Thus, in equation (2.33), only the first term in the sum need be considered. The 
electron number density is given by 



P 



° 23/27r2A3 Pq VmocV ^'''^\kTj 
where we have kept only the 1/ — term. For p < pB, the electron Fermi temperature is 



Tp = Ev/k = 2.70 Pil, (FcPi)' K, 



(for p < pb) 



(2.37) 
(2.38) 
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where pi = p/(l g cm~ ). Because electrons are confined by the magnetic fields to smaller volumes, electron degeneracy occurs 
at higher densities than in the absence of the magnetic field. Note that for T ^ Tp, equation (2.33) reduces to Pc ~ n^kT 
even when the field is strongly quantizing (p < p_g). The ion pressure, on the other hand, is well approximated by 

P.„^^. (2.39) 

The Coulomb correction is given by (e.g., Shapiro & Teukolsky 1983) 

47r\i/3 



Pc. 



loul 



-l^^fyZ^WJ^. (2.40) 

The equation of st ate, i ncluding the effects of degeneracy, is determined by substituting the different components of the 
pressure into equation ( |2.32| ) to obtain a relation between a given total pressure P{t) and temper ature T{t). The root of this 
equation yields the degeneracy parameter r;(r) — pc/(fcT) which is then substituted into equation ( 2.37 ) to obtain the electron 
number density ric and the total density p = {Amp/Z)nc. We use the rational function approximation given by Antia (1993) to 
evaluate the Fermi-Dirac integrals. For all the nonmagnetic and magnetic atmosphere models considered here, the difference 
between the densities calculated from the ideal gas equation of state and the equation of state including electron degeneracy 
is only a few percent in the innermost layers, and the effect on the surface spectrum is negligible (see Section 



2.5 Opacities 



In a magnetized electron-ion plasma, the scattering and free-free absorption opacities depend on the direction of propagation 
k (more precisely, on the angle 9b between k and 13) and the normal modes {j = 1,2) of the electromagnetic waves. In 
rotating coordinates, with the z-axis along B, the components of the polarization vector e-* are given by (see Appendix^) 



j 1 2 
J |2 



V2 



2(1 + Kf+Kl^, 
{Kj sin^B — Kz,j cos 6b)^ 



1 ± {Kj cos 9b + Kzj sin t 



where 



1/2 



UeV (^1 — Ui — M ^) sin6'_B cos^sifj 



1/2 . , 

' Me V Sm f 



p = 



(1 - itc) (1 -ui) -V [(1 - ui) (1 - Mo cos2 9b) - Mui sin^ ( 

1 + M^ 



1/2 • 2/1 
Me sm Ob 

2(1 - v) COS 9b 



M 



and M = Ara^j (Znic). In equations (2.4^ )-( p.45| ), we have defined 

.2 , .2 



Mc 



2 



(2.41) 
(2.42) 

(2.43) 

(2.44) 
(2.45) 

(2.46) 



the electron cyclotron frequency cjse, the ion cyclotron frequency ojsi, and the electron plasma frequency ix)p are given by 
eB 

fiujB^ = = 11.58Bi2keV 

rrieC 

huBi ~ ft = 6.305 -Bi2 ^ ^ 



(!) 



eV 



hujp = 



47re ric 



1/2 



V\ 1/2 
Z\ ' 1/2 



= 28.71 ( - ) pi'^ eV, 



(2.47) 
(2.48) 

(2.49) 



respectively. Note that in the coordinate frame where the z-axis is along k and the a;-axis is in the k-B plane, the polarization 
vector can be written as (see Appendix ^) 



-.{iKj,l,iK.,j). 



(2.50) 



Thus the modes are elliptically polarized with Kj being the ratio of the axes of the ellipse. It is easy to show that Ki = —K2^ 
and \Ki\ < 1; this implies that the extraordinary mode {j — 1) is mostly perpendicular to the k-B plane, while the ordinary 
mode {j — 2) lies mainly in the k-B plane. For lu <^ ujbc, the modes are nearly linearly polarized at most angles except close 
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to parallel propagation 6b = Q and vr, where the polarization is circular since Kj 



-ly . In this limit, 



rJ 1 2 



and 



and le^l ~ (3 [see Section 2.6, specifically eqs. (2.64) and (2.65) 



It is evident from equation (2.44) that the component of e-' along k is of order v oc p/ui^, and thus, at sufficiently 
low densities, Kz,j can be neglected so that the modes are transverse. Previous studies adopted this "transverse-mode" 
approximation when treating the radiative opacities in a magnetized medium (e.g., Pavlov et al. 1995). As we shall see in 
Section ^, since the radiation can directly emerge from deep layers where « ^ 1, the transverse-mode approximation cannot 
always be justified. 

In a magnetized medium, the electron scattering opacity from mode j into mode i is given by (Ventura 1979; Ventura, 
Nagel, & Meszaros 1979) 



1 



p ^ (OJ + QLJfle)^ + vl 
a — — 1 



J |2 



4' 



(2.51) 



where is the angle integral given by 



4* — — 



dk lei 



The electron scattering opacity from mode j is 

1 2 
OS '^cO'T I 7 |2 . 



a — — 1 



(2.52) 



(2.53) 



where Aa = Y^\ 
relation 



J. In the transverse-mode approximation (Kz,j = 0), the polarization vector e-* satisfies the completeness 
= (l +cos^S_b) /2 and ^ jejp = sin^^s, and thus Aa ~ 1. Note that a resonance develops near 
the electron cyclotron frequency for the circularly polarized mode with its electric field rotating in the same direction as the 

■ ^0 



4P 



electron gyration; however, this resonance appears only in the extraordinary mode since one can easily show that |el 
at = LOBe- For the photon energies of interest in this paper, i.e., lo <^ lobc, a suppression factor (ujjwBe)^ in the opacity 
results from the strong confinement of electrons perpendicular to the magnetic field. Similar features appear in the electron 
free-free absorption opacity (e.g., Virtamo & Jauho 1975; Pavlov & Panov 1976; Nagel & Ventura 1983) 



P 



1 

V 

a--\ 



I i i2-ff 



(2.54) 



where 



Qo 



- 2^2 3 ft^c^ (2m^\^l'^ n^m 



) =ao 



(2.55) 



ap = /{he) is the fine structure constant, and OfQ and (f^ are the free-free absorption coefficient and velocity-averaged 
free-free Gaunt factor, respectively, in the absence of a magnetic field. In equation (2.54), — and go = g^ are the 
velocity-averaged free-free Gaunt factors in a magnetic field, which we evaluate using the expressions given in Meszaros (1992) 
(see also Nagel 1980; note these expressions for the Gaunt factors neglect transitions to the excited electron Landau levels; 
see Pavlov & Panov 1976 for the complete expressions). 

In an ionized medium, the ions also contribute to the total scattering and absorption opacities. Analogous to equa- 



tions (2.51)-(2.54), we have 



Z rric 
Arriry 



Arriry 



riiCTT 



1 

a— — 1 
1 

a = -l 



Q-LUBi)^ + V? 



aujBi)^ + v'f 



\e-L\ Ao 



(2.56) 



(2.57) 



fF,i 



Arrir, 



{lu — aLUBiY + v\ 



leal goc- 



(2.58) 



Note that the ion cyclotron resonance occurs for a = -f 1, i.e., when the electric field of the mode rotates in the same direction 
as the ion gyration. The total scattering and absorption opacities in a fully ionized medium is then the sum of the electron 
and ion component s, namely ti ^^ = k^^ -\- and k^'^" = 



In equations ( 2.51 )-( 2.58 ), we have included damping through 



(2e /SiricC and Vr.i = {Z mc/Amp)i^r,c are radiative damping rates and i^c.o 



h !^c,c and i/i — i/rj + i/c 
(aog*/nocrT)j^r,c and i/c,i 



where fr.e ~ 
{■mc/Amp)uc,c 



are coUisional damping rates (see Pavlov et al. 1995 and references therein) . Note that the natural width of the ion cyclotron 
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resonance is determined by the dimensionless damping rate 71 = Vi/uj — 7r,i + 70,4, with radiative and colhsional contributions 
given by 



7r,i 



5.18 X 10" 



~A VkeVy 



7,,,; = 8.42 X lO-'^-^sS 



l(f_K 

T 



Hi 



/keVy 
lO^"* cm-^y \huj J 



Also, the thermal width of the ion cyclotron line width is of order 
/ 2kT 



7th 



\ ArUpC^ 



1.46 X W'-'A 



VkeV; 



(2.59) 
(2.60) 

(2.61) 



the Doppler broadening would presumably change the line profile from the Lorentzian to the Voigt profile (e.g., Mihalas 1978). 
Since 7r,i, 70,1, 7th ^ 1, damping is negligible except very near resonance (see Section 2.6). 



2.6 Behavior of Opacities and Ion Cyclotron Resonance 

Figure ^ shows the angle- averaged absorption opacity Kf°^ and scattering opacity K^'^ = X/i=i I®^® ^1®- ( 2.25)-( 2.27)] 
as functions of energy for magnetic fields B — 4.7 x 10^^ G and lO'^* G at densities and temperatures characteristic of NS 
atmospheres. For angles 6b not too close to 0° or 180° (e.g., 10° ^ 9b ^ 170°), the unaveraged opacities and i^f^^ resemble 
KY and Kf°^ in qualitative behavior, while exhibiting strong angle dependence for 6b near 0° or 180°. The most prominent 
feature in Fig. ^ is the ion cyclotron resonance at TiLOBi ~ 0.63 -B14 keV (for Z = A — 1), where Bn = _B/(10^^ G), for the 
extraordinary mode. Note that the cyclotron feature does not possess the usual Voigt profile that characterizes spectral lines; 
the feature is rather broad with AE/E ~ 1 + 10"'^ B^^pi [see eqs. (2.69)-(p.70|) below], which is far wider than the natural 



or thermal width [see eqs. (2.5t )-(2.61)]. Away from the resonance, the extraordinary-mode opacities are suppressed relative 
to the zero-field values by a factor (ujjbJBe)^ (see below). For the ordinary mode, the opacities are similar to the nonmagnetic 
values (see below) except at high densities/low frequencies, i.e., « <; 1, where the effects of A'zj become important [see 
eq. 

The behavior of the opacities, especially the ion cyclotron feature, shown in Fig. |l| can be understood as follows. For the 
magnetic field strength and photon energy of interest in this paper, luo <C /loiBe = 11.6 -B12 keV, the scattering and free-free 
absorption opacities can be written as 



p 



(2.62) 



with 



.J 1 2 



+ 



M''- 



(1 + V^)' 



(2.63) 



where M' = {Amp/Zmc)Z for scattering and M' = {Amp/Zma)Z^^'^ for absorption, and we have set the free-free Gaunt 
factors to unity in n'f"^ and Ac = 1 in for simplicity. In the following, we shall consider photon energies such that 



J 

\u\ — Ij <; 7i, 7th [see eqs. (2.59)-(2.61)], and we shall drop the 7? term from equation (2.63). Also for simplicity, we shall set 
Z = A = 1. 

For 6b not too small, i.e., sin'^ 6b/ cos 6s ~ 1, the polarization parameter (3 [see eq. ( 2.45 )] satisfies 3> 1 except when 
u) is very close to the ion cyclotron resonance (i.e., \1 — u\ — {1 + v)/M\ ^ u^^^^) or when hu) ^ hUp/LUBe = 0.07^15^""^ eV. 
In the fimit of |/?| > 1, we find Ki = -1/(2/3) and K2 = 2/3, and thus 



= - 1 =F — cos ( 



,2,2 1 2 , 
le±l = 2 '^^^ ' 



1 ± 



/3 cos t 



1 



1 



sin 6b, 



where, for simplicity, we have taken the transverse approximation so that K^,] 
we obtain, for |/3| 3> 1, 



(2.64) 
(2.65) 



0. Substituting these into equation (2.63), 



1 + 



(1 



■ cos^ 6b/ sm^ 6b + Ui{l 



1 



cos 6b + Ui sin 



- Ui) 



2 Ui{l + u{) 

cos 6b — ^ 



(1 



-vf/sm^ 



(2.66) 
(2.67) 



Clearly, away from ion cyclotron resonance, 



OC lie ^ 



(bj / LOBiif" , and thus the opacities for the extraordinary mode are 
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suppressed compared to the zero-field values: 



Ki OC 



(- 



abs 

OC 



1 



Be 



~riLj/kT\ 



(2.68) 



At high densities/low frequencies, such that w <; 1, these are multiplied by a term of order 1 + 0{v'^); for u ^ 1, we have 
Ki^ OC u}p/{uj%^Lo^). For the ordinary mode, ~ 1, and thus the opacities are approximately unchanged from the zero-field 
values. 

Equations (2.66) and (2.67) also apply to photon energies near the ion cyclotron resonance as long as |1 — Ui| ^ {l + v)/M, 



so that |/3| ^ 1 (note that we are considering typical angles where | cos^_gj ~ 1). It is evident from equation (2.66) that, for 
the extraordinary mode, the opacities exhibit a broad peak/feature around uisi, i.e., the term contained inside the brackets 
in equation (2.66) becomes significantly greater than unity. The width of the feature, |1 — Uij, is of order 

1 1/2 



where 



UJ — LOBi 




OJBi 













2 + (l 



2 cos 



sm 



B 



.4.6 X 1012 G 



(2.69) 



(2.70) 



It is important to note that, although the natural and thermal width of the ion cyclotron line are rather narrow [see eqs. ( |2.59| )- 
(2.61)], the cyclotron feature as described by equations (2.66) and (2.69) is broad, i.e., the opacities are significantly affected 
by the ion cyclotron resonance for 0.3 ^ uj/uJBi ^ 3. Also, equation (2.67) indicates that, although the ion cyclotron resonance 



formally occurs for the ordinary mode, its strength is diminished by the {uj/ujbr) = {me/m-p) suppression. 

Finally, we note that the presence of ions induces a mode-collapse point very near uobi- This can be seen from equa- 
tion (2.45), which shows that /9 = when 1 — Wi — (1 -I- v)/M = or w = uJBi[l -I- (1 + v)/2M]. At this point, the two photon 
modes become degenerate (both are circularly polarized), and the transport equation (2.1) formally breaks down. This ion- 
induced mode collapse is independent of the well-known mode collapse due to vacuum polarization (Pavlov & Shibanov 1979; 
Ventura et al. 1979; see Section]^. Note, however, that the ion-induced mode-collapse feature is very narrow [|/3| ^ 1 requires 
|1 — — (1 + v)/M\ -C Ue^^'^], so it is completely "buried" by the much wider and more prominent ion cyclotron feature. 



3 NUMERICAL METHOD 

3.1 Magnetic Atmospheres with Full Radiative Transport 

When presenting results of atmosphere models based on the full transport equations, we consider only the case of B per- 
pendicular to the surface (Qb = 0). Thus the radiation intensity II depends on depth r, energy E = Tiuj = hv, and angle 

e. 

We construct a grid in Thomson depth, energy, and angle with intervals that are equally spaced logarithmically from 
10~* to ~ 10^ in depth and 10"^ to ~ 10 keV in energy and spaced every five degrees in angle. For typical calculations, six 
grid points are used per decade in depth, twelve grid points are used per decade in energy, and nineteen total grid points in 
angle. We also increase the energy resolution near the ion cyclotron resonance. 

The temperature profile is initially assumed to be the grey profile, 

1/4 



T(rR) = Toff 



The Rosseland mean depth tr is given by dra 
1 37r 1 , dB^ 



/ Ko^)dT, where the Rosseland mean opacity is 



4crsBr^ 



dT 



■du 



(3.1) 



(3.2) 



and Ij is given by equation (2.23). Since Ij and k are functions of p and T themselves, the profile (3.1) must be constructed 



by an iterative process: the temperature profile is first taken to be T{t) = refr[3/4(r + 2/i)Y^'^; 
hence and kP', are then determined, which allows the calculation of tb.{t)\ using equation (3.1) with this tr, a revised 



the density profile, and 



temperature profile is obtained; the process is repeated until the new ni{T) is not significantly different from the previous 
TR,(r). To achieve |Ar/r| and |Atr,/tr,| < 0.1%, three iterations are required. 

To construct self-consistent atmosphere models requires successive global iterations, where the temperature profile T{t) 
is adjusted from the previous iteration to satisfy radiativ e equ il ibrium . At each global iteration, the profiles T{t) and p(r) are 
consid ered fixed. We solve the RTE ( ^.8| ) [see also eqs. ( 2.15 )-( 2.17 )], together with the boundary conditions [see eqs. (^.Sj)- 
(2.1C)] by the finite difference scheme described in Mihalas (1978). This results in a set of equations which have the form of 
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a tridiagonal matrix that can be solved by the Feautrier procedure of forward-ehmination and back-substitution. We have 
also implemented an improved Feautrier method to reduce errors arising from machine precision as suggested by Rybicki & 
Hummer (1991). Note that, in the presence of scattering, the source function 5*^ [which involves mode coupling; see eq. (2.16)] 
is not known prior to solving the RTE. We adopt the strategy that, for each global iteration, SI is calculated using the 
solution for ii from the previous global iteration, with Si — Bu/2 for the first global iteration. The full RTE is solved with 
this source function to obtain ii, and the specific flux is calculated from finite differencing equation (2.7). Note that, unless 
the temperature correction from the previous iteration is relatively small, so that the source function calculated from the 
previous ij approximates the true current source function, the global iteration may not converge.^ 

The current solutio n does n ot, in general, satisfy the radiative equilibrium constraint nor does it yield a constant flux 
at every depth [see eq. (2.13) and ( [2.14| ), respectively]. To satisfy both conditions requires correcting the temperature profile. 
We use a variation of the Unsold-Lucy temperature correction method as described in Mihalas (1978) but modified to account 
for full radiative transfer by two propagation modes in a magnetic medium. The radiative equilibrium constraint given by 
equation (|2.13|) can be written as 



dr 



ti p 

u 

i-CS 

^0 



(3.3) 



where u = j duui is the total radiation energy density, if — aT^ = (4(7sb/c)T'' is the blackbody energy density, and 

K'' and are the absorption mean opacity and Planck mean opacity defined by 



-E 



dv 



2 



j = l 



dk «:f=(k)ii(k) 



dv \ dkKf"(k)ttJ; 



(3.4) 



(3.5) 



respectively. On the other hand, equation (2.7) can be integrated out formally to give ii in terms of jl, and we then have 



cu(r)=2y \ dv \ dkii(r)= / dr'^^^^F,(r') + 2F,(0), 
—"^J J+ Jo 



(3.6) 



where we have used cu{0) 
defined by 



2Fz(0), Fz{t) is the total flux at depth r [see eq. (2.14)], and k, (t) is the flux mean opacity 



2 



E 



du 



dk-4°*(k)/^(k). 



(3.7) 



Combining equations (3.3) and (3.6) gives an expression for 'uF\T{t)\ in terms of Fzir). The desired temperature correction 
Ar(r) needed to satisfy radiative equilibrium is then 



Ar(r) 



160-SBr(T)3 

where AF^r) = as^T^s 
Ar(r) ^ 



.P(r) 



d[AF.(r)] AcJ(r) 
AcP(r) 



dr 



dr 



AF,{r') + 2AF,{0) 



16o-sBr(r)3 1 «:P(r) 



Fz{t). Using equation (3.3) to replace dF^/dr in equation (3.8), we obtain 



[«-'(r)u(r)-KP(r)uP(r)] + 



"dr'^AF.(r') 



2AF40) 



(3.8) 



(3.9) 



), corresponds to the temperature correction in the A- 



Note that the flrst term in equation (3.9), which is oc 
iteration procedure ATa = (k^u— K,^u^){du^ /dT)~^ and is important near the surface but is small in the deeper layers where 
the atmosphere approaches a blackbody, while the remaining terms provide corrections in these deeper layers. In practice, 
we find that the Unsold-Lucy temperature correction method as defined by equation (3.9) tends to over-correct; therefore 
we only use ~70-90% of the temperature correction. The process of determining the radiation intensity from the RTE for a 
given temperature profile, estimating and applying the temperature correction, and then recalculating the radiation intensity 



i An alternative is to calculate Sl iteratively within each global iteration. The source function during each global iteration is first 
taken to be Si = By 12, and thus the two modes are initially decoupled. The RTE is then solved for each mode j. The resulting ii is 
used to determine a revised source function, which can then be reinserted into the RTE. This procedure is repeated until the the source 
function converges, e.g., \C\SH Si\ < 5%. Obviously, this method requires longer computation time. In addition, we have implemented a 
procedure that reformulates the RTE using Eddington factors (see, e.g., Rybicki & Lightman 1979 for the analogous nonmagnetic RTE 
using Eddington factors), but this does not lead to significant reduction in computation time. 
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is repeated until convergence of the solution is achieved. We note here that, to decrease the number of iterations required 
for convergence and the total computation time, we can use as the initial temperature profile the result obtained from the 
diffusion calculation (see Section 3.2); such an initial temperature profile would already be close to radiative equilibrium. 

Three criteria are used to indicate convergence. The first is that the temperature correction between successive global 
iterations becomes small, e.g., |AT/r| ^ 1%. The second is to check that radiative equilibrium [see eq. (2.f3)] is sufficiently 
satisfied, or {hC^u— v^) /{k^ u) is sufficiently small [see eq. (3.3)]. The third criterion is that Fz is sufficiently close to ctsbT^j 
[see eq. (pT^J. Figure I shows an example of these convergence tests after about twenty global iterations for our full radiative 
transport model of a fully ionized, pure hydrogen atmosphere with Tcff = 5 x fO'' K, B — fO^* G, and Os = 0. Further 
iterations reduce the deviations, though the convergence is slower since the temperature corrections are already small. The 
numerical results presented in Section ^ have all reached the convergence level similar to Fig. |^. 



3.2 Magnetic Atmospheres with Diffusion Approximation 



Using the diffusion approximation of radiative transport (see Section 2^), we construct magnetic atm osph ere models for 



general Qb- We fo llow a similar procedure as the full radiative trans fer m ethod desc ribed in Section 3.i. We solve the 
approximate RTE ( 2.24 ) with boundary conditions given by equations ( 2.2^ ) and ( ^.29| ) using the Feautrier method, which 
involves only depth and photon energy grid points. The modified Unsold-Lucy temperature correction can be derived using 
the specific flux perpendicular to the stellar surface, F^ .^ ~ -■clj{du-l/dz), and equation (2.30). The resulting AT(r) takes 
the same form as in equations (3.^) and (3.!-), except that the mean opacities k'', , and kF are given by 



2 

-E 



2 

—y 

2 

-y 

F. ^ 



dr^abs 7 



T r^abs P 

au Kj 



dv—Fi, 



(3.f0) 



(3.1f) 



(3.12) 



respectively. The expressions for k'' and are essentially the same as in Section 3.1. Atmosphere models for B < 10^"^ G 



based on the diffusion approximation have been constructed by Shibanov et al. (1992); Pavlov et al. (1995) using a different 
iteration scheme, i.e., A-iteration. 



3.3 Nonmagnetic Atmospheres 

For comparison, we also construct nonmagnetic atmosphere models by solving the standard RTE (both full transport and in 
the diffusion approximation) using a similar method as described in Section 3.1 and Section 3.2, In this case, the opacities are 
just due to nonmagnetic electron scattering and free-free absorption. We implement the RTE with Eddington factors to more 
accurately determine the source function. We use the standard Feautrier method (with improvements described in Rybicki & 
Hummer 1991) and the standard Unsold-Lucy temperature correction scheme described in Mihalas (1978). 



4 NUMERICAL RESULTS 
4.1 Atmosphere Structure 

Figure ^ shows the temperature proflle as a function of Thomson depth r and Rosseland depth tr for hydrogen atmospheres 
with Tcfi = 5 X 10^ K, magnetic fields B = 4.7 x 10^^ G and lO^'' G, and the magnetic field oriented perpendicular to the 
surface (Gg =0). For comparison. Fig. |^ also shows the temperature profile of our nonmagnetic atmosphere model and the 
grey profile given by equation ( |3.l[ ). Note that, since the Rosseland mean opacity [eq. (3.2)] is dominated by the (small) 



extraordinary-mode opacity, the Rosseland depth tr{t) is significantly smaller at high fields than at _B = 0. It is evident that 
all models approach the grey profile a.t tr ^ 1, although this occurs deeper in the atmosphere for models with higher magnetic 
fields. At small depths of r 1, the gas temperature is smaller, by a factor of a few, than predicted by the grey profile. This is 
because, for r ^ 1, the extraordinary-mode photons have already decoupled from the matter, and the temperature is primarily 
determined by the ordinary-mode absorption opacity, which is strongly nongrey oc approximately). The plateau in 

the r(r) profile at r ~ 1 - 100 in the B = lO" G model arises because the extraordinary-mode photons decouple from the 
matter at r ^ 100 (see Section ^^); such a plateau is more distinctive for high magnetic fields and becomes less pronounced 
for smaller magnetic fields. 
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Figure ^ shows the temperature versus density profiles for various atmosphere models considered in this paper (see 
Section 4.3 for the spectral properties of these models). The plateau in temperature for the B — 10^* models correspond to 
the similar feature shown in Fig. ^ Note that, in all models, the electrons are nondegenerate at the relevant densities, i.e., 
T ^ Tf [see eq. ( ^.38| )]. Thus, even for high (quantizing) magnetic fields, the equation of state deviates from the classical 
ideal gas law by less than a few percent (see Section 2.4). At a given B, the temperature profile, T{t), of the He atmosphere is 
very close to that of the H atmosphere (see Fig. this is expected since, apart from the ion cyclotron resonances, the H and 
He opacities have the same frequency dependence, and differ only by a factor of ~ 2 in magnitude. With P ~ (1 + Z)pkT/A, 
we see that the He atmosphere has a higher density, by a factor of ~ 16/3, than the H atmosphere at the same temperature 
(see Fig. a. 



4.2 Photon-Matter Decoupling Region 



To better understand the atmosphere structure and the emergent spectra (see Section 4.3), it is useful to develop a rough 
idea of the physical conditions in the region of the atmosphere where photons of different modes and energies decouple from 
the matter. Because of the angle dependence of the opacities, the exact location of the "decoupling layer" is ambiguous. 
One possible approach is to use the averaged opacities defined in Section 2.3; the effective opacity for photon-matter energy 
exchange is Kf = {Kj^f-f^'^ (see Rybicki & Lightman 1979), where Kf^ = {plj)'^ and Ij and Kf'' are given by 



equations (2.23) and (2.25), respectively. In fact, if mode coupling due to scattering is neglected, equation (2.24) can be 
approximately integrated out to yield ul{y = 0) cx u^{y ~ l/Kf^), i.e., the emergent radiation density (of mode j) is 
determined by the blackbody radiation density evaluated at the depth where the column density y ~ 1/ Kf^ . We therefore 
define the effective optical depth 



dr' 



tot -J 1/2 ^ 



dr' 



K 



abs \ 1/2 



(4.1) 



so that photons of mode j and frequency p decouple from the matter at ~ 1. Figure ^ shows the Thomson depth 
and local temperature and density as functions of photon energy at the decoupling layer for the H atmosphere model with 
Tcff = 5 X lO*' K, Qb = 0, and B = 4.7 x 10^^ G and lO" G (see also Fig. 1). Note that the efltective depth defined in 
equation ( |4.l[ ) errs for photons propagating along the magnetic field since the opacities are generally lower near 9b = 0; 
hence, these photons decouple from the matter at deeper layers than those indicated in Fig. |^. Nevertheless, Fig. ^ gives the 
typical conditions where the observed photons of energy E are generated. It is clear that the extraordinary-mode photons 
can emerge from deep in the atmosphere where plasma effects on the opacities are not negligible (see Section |2.5| and Fig. |l|) . 
Photons with energies near the ion cyclotron resonance EBi ~ 0.63 B14 keV decouple in the lower temperature region; this 
gives rise to the absorption feature in the emergent radiation spectrum (see Section 4.3). 



4.3 Spectra 

Figure ^ shows the spectra of the hydrogen atmospheres with TcS — 5x lO" K, Os =0, and magnetic fields of B = 4.7 x 10^^, 
10^*, and 5 x 10^* G. Also plotted are the spectrum of our nonmagnetic fully ionized H atmosphere model and the blackbody 
spectrum at T = 5 x 10^ K. The spectra clearly show a significantly harder high energy tail and a depletion of low energy 
photons relative to the blackbody. For example, in the case of _B = 10^* G, the fiux at iJ = 5 keV is an order of magnitude 
larger than the corresponding blackbody flux. At a given energy E (and for the same Tcff), the flux depends on the field 
strength in a non-monotonic fashion for the models depicted in Fig. ^. For example, at iJ = 7 keV, the flux increases in the 
order of B = 4.7 x 10^^ 10^*, 0, 5 x lO" G. The hard tails were already noted in previous studies of NS atmospheres with 
B ^ 10" G (e.g., Pavlov et al. 1995; see also references cited in Section 111) , and they arise because high energy photons have 
smaller opacities and thus decouple from deeper, hotter layers (see Fig. It is evident that fitting the high energy tail of 
observations with a blackbody will yield overestimates of the surface temperature of the NS. 

Most prominent in the spectra of highly magnetized atmospheres is the absorption feature at the ion cyclotron resonance 
Ebi ~ 0.63 [Z/A)B\A keV. This can be directly traced back to the corresponding resonance feature in the opacities (see 
Fig. ^). The ion cyclotron feature is broad, with AE/Ebi ~ 1, and is much broader than the natural or thermal width of the 
line [see Section 2.6 for a discussion and eq. (2.6S;) in particular]^. We find that the ion cyclotron feature is most pronounced 



when it appears beyond the blackbody peak, i.e., when Ebx <; SfcTcft . For the models depicted in Fig. ^ the equivalent width 



S In discussing the cyclotron feature in this paper, we are referring to the broad depletion of flux around the ion cyclotron energy rather 
than the yp.rv narrow {HE/E ~ 10~^) natural/thermal line; such a line has negligible equivalent width and is not always resolved in 

Figs. 0-1 In} 
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of the Esi = 3.15 keV line at Bn = 5 is about 1.9 keV, the 0.63 keV line at B14 = 1 is 90 eV, and the 29.6 eV line (not 
shown) at B12 ~ 4.7 is only 1.3 eV. 

Figure ^ shows the spectra of H atmospheres at two different effective temperatures: Toff = 10^ K and 7 x 10^ K. The 
spectra are similar to those of the Tcs = 5 x 10® K models shown in Fig. ^ The cyclotron resonance at 0.63 keV (for B14 — 1) 
becomes prominent for TcH — 10^ K since the condition Ebi 3 kTcB is satisfied. Note that the TcS ~ K models should 
be considered for illustrative purposes only since we expect the opacities to be significantly modified by neutral atoms and 
even molecules at such a low temperature (see Section ^ . 

As noted before, in strong magnetic fields, the thermal radiation is primarily carried by the extraordinary-mode photons. 
This can be seen clearly in Figure |^. Obviously the emergent radiation is expected to be polarized. Also note that the ion 
cyclotro n re sonance feature occurs only for the extraordinary mode, as expected from the behavior of the opacities (see 
Sections |2.5| a nd ^ ) . 

Figure H shows the spectra of the He atmospheres with T^g = 5 x 10® K, Os = 0, and magnetic fields of _B = 0, 4.7 x 10^^, 
10^* and 5 x 10^* G. These should be compared to the spectra of H atmospheres shown in Fig. ^ Although the He atmosphere 
is denser than the H atmosphere (see Section and Fig. there is no clear distinction between their spectra except for 
the ion cyclotron resonance features. For example, at B = 10^* G and T^s — 5 x 10® K, the fiux of the He atmosphere at 
E = 6 keV is lower than that of the H atmosphere by only a factor of 1.2. Clearly, it would be difficult to distinguish a 
H atmosphere from a He atmosphere based on spectral information in the 1 — 10 keV band alone. We also note that the 
equivalent widths of the ion cyclotron features for the helium atmospheres are about 0.51 keV for the Ebz ~ 1.58 keV line at 
Bi4 = 5 and about 10 eV for the Ebz ~ 0.32 keV line at B14 — 1. 



4.4 Full Radiative Transport Versus Diffusion Approximation 



The atmosphere models presented in Sections 4.1-4.3 are based on solutions of the full radiative transfer equations (see 



Sections 2.1 and 3.1). These models are currently restricted to the case where the magnetic field is perpendicular to the 



stellar surface, i.e., Ob ~ 0. Using the diffusion approximation (see Sections |2.3| and 3.2), we construct atmosphere models 
for general Ob's. 

Figure |lo| shows the temperature-density profiles of the H atmosphere models with B = lO^"* G, Tofi = 5 x 10® K, and 
Ob = and 90° based on the diffusion approximation. These are compared with the result of the corresponding full transport 
model with Ob = 0. Clearly, in the diffusion models, the surface temperature (where p ^ 0.01 g/cm'') for Os = 90° is 
significantly less than for 0b ~ (by roughly a factor of 1.4 for the models considered here); this was already noted by 
Shibanov et al. (1992) and Pavlov et al. (1995) for the B ~ 10^^ G models. At lO'^* G, the temperature plateaus (due to the 
decoupling of the extraordinary-mode photons; see Section are also different for 0s = and Ob = 90°. The difference 
in the temperature profiles for different magnetic field geometries implies that different parts of the NS surface will have 
different temperatures, e.g., the magnetic equator is cooler than the magnetic poles, even if T^b is uniform. Clearly, a full 
understanding of the NS surface temperature distribution requires detailed modeling of the three-dimensional atmosphere 
and of the anisotropic heat transport through the NS crust; this is beyond the scope of this paper. Comparing the diffusion 
model with the full transport model, we see from Figure ^that, at Os = 0, the temperature in the full model is significantly 
lower than that of the diffusion model in the outermost layers of the atmosphere. 

Figure [ll| shows that, despite the difference in the temperature profiles, the spectra of the diffusion models for Os = and 
90° and the spectra of the full transport model are rather similar. For B14 = 1, the equivalent width of the cyclotron feature 
at Esi = 0.63 keV is about 80 eV for ©3 = and 50 eV for Ob = 90°, while for Bu = 5, the feature at Egi = 3.15 keV 
has a width of about 1.9 keV for ©b = and 90°. We also find that the polarized emission is stronger for Ob ~ 90° than for 
Ob = 0, as already noted by Shibanov et al. (1992) for the B = 4.7 x 10^^ G models. 



5 DISCUSSION 

In this paper, we have constructed models of magnetized neutron star atmospheres composed of ionized hydrogen or helium. 
We focused on the superstrong field regime [B > lO^** G) as directly relevant to magnetars (see Section |^). We solve the full, 
angle-dependent radiative transfer equations for the coupled polarization modes and include the ion cyclotron and plasma 
effects on the opacities and polarization modes. Since the magnetic field greatly suppresses the opacities of the extraordinary- 
mode photons, the thermal radiation emerges from deep in the atmosphere (p ~ lO'^ g/cm"^ for B — lO^** G), where the 
medium effect is non-negligible. As already noted in previous works on neutron star atmospheres with moderate (~ lO^'^ G) 
or negligible magnetic fields, the thermal emission from a highly magnetized neutron star is polarized, with the spectrum 
harder than the blackbody. We find that the ion cyclotron resonance, at EB^ = 0.63 {Z/A){B/W^* G) keV, can give rise to a 
broad (AE/E ~ 1) and potentially observable absorption feature in the spectrum, especially when i?Bi ^ 3 kT^s- Clearly the 
detection of such a feature would provide an important diagnostic for magnetars. 
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We note that the spectra presented in this paper (see Figs. and |ll|) correspond to emission from a local patch 
of the neutron star surface. To directly compare with observations, one must calculate synthetic spectra from the whole 
stellar surface, taking into account the effect of gravitational redshift and light-bending. To do this, one must know the 
distribution of the magnetic field (both magnitude and direction) and the effective temperature over the stellar surface. For 
a given field geometry, the surface temperature distribution may be obtained from crustal heat conduction calculations (e.g., 
Hernquist 1985; Schaaf 1990; Page 1995; Heyl & Hernquist 1998), but this is intrinsically coupled to the properties of the 
atmosphere; the problem is further complicated by the existence of lateral radiative flux (see Section ^.3[) and the nonuniform 



temperature distribution induced by the magnetic field in the atmosphere (see Section 4.4). Obviously we expect the cyclotron 
feature to be broadened and less deep if different parts of the surface with different field strengths contribute similarly to the 
observed flux (see Zane et al. 2001 for a calculation based on a dipole field and an approximate temperature distribution). 

Several physical effects have been neglected in our atmosphere models (see Ozel 2001; Zane et al. 2001, for related 
works on magnetar atmospheres; see also Section ^ for references of previous works). We have assumed that the atmosphere 
is fully ionized. The problem of ionization equilibrium in a highly magnetized medium is a complicated one owing to the 
nontrivial coupling between the center-of-mass motion and the internal atomic structure and to the relatively high densities 
of the atmosphere (see Lai 2001 and references therein). While estimates based on the calculations of Potekhin et al. (1999) 
indicate that the fraction of neutral atoms is no more than a few percent for T ^ 5 x lO'' K at B ~ 10^* G, it should 
be noted that even a small neutral fraction can potentially affect the radiative opacities. For example, the ionization edge 
of a stationary H atom (in the ground state) occurs at energy Ei ~ 4.4 (In 6)^ eV, where b = 5/(2.35 x 10^ G) (thus, 
Ei = 0.16, 0.31, 0.54, 0.87 keV for B = 10^^ lQl^ lO^"*, 10^^ G, respectively; see Lai 2001). The ratio of the free-free 
and bound-free opacities (for the extraordinary mode) at E = Ei is of order 10~^' pi/{Tg^^ Bi4fH), where pi — p/{l g cm^''), 
Te = r/(10'' K), Bi4 = B/(10^* G), and /h is the fraction of H atoms in the "centered" ground state (see Lai 2001). Thus 
one could in principle expect atomic features, broadened by the "motional Stark effect," in the spectra; these could blend 
with the ion cyclotron feature, at Ebz, studied in this paper. 

Another caveat of our models is the neglect of the vacuum polarization effect. In a strong magnetic fleld, vacuum 
polarization (in which photons are temporarily converted into electron-positron pairs) contributes to the refractive index by 
a term of order S^p = {ap /457v){B / Bq)^ , where ap = 1/137 and Bq = 4.4 x 10^^ G, and the effect tends to induce linear 
polarization of the modes (Adler 1971). Thus for 35vp ^ (tOp/uj)^ (where tOp is the electron plasma frequency), or, for photon 
energies hu) ^ TuVvp = (157r/aF)^/'^(tjp/tJse)mcC^ = 1.0 p\^^ B^^ keV^], the normal modes (and thus the radiative opacities) 
are modifled by vacuum polarization (e.g., Gnedin et al. 1978, Meszaros & Ventura 1979). Previous works (see Meszaros 1992 
for a review; see also Shibanov et al. 1992, Ozel 2001 for atmosphere models which include the vacuum polarization effect; 
note that Ozel 2001 studied the B 3> -Bq regime but used the vacuum polarization expressions valid only for B ^ Bq) 
neglected the response of the ions in calculating the dielectric tensor. However, as shown in this paper, for B ^ 10^* G, the 
ion cyclotron resonance gives rise to a prominent feature in the polarization mode vector and in the opacities. Therefore one 
must consider the combined system of electrons, ions, and polarized vacuum in order to model the spectra in the superstrong 
field regime (Ho & Lai 2001). Also, vacuum polarization causes the opacities to exhibit a resonance feature ai E = huj^p 
(Pavlov & Shibanov 1979; Ventura et al. 1979; see Meszaros 1992 for a review). Since huvp cx p^^^ and since atmosphere 
models necessarily require discrete depth (density) and energy grid points, such features can appear in the model spectra, 
while, in a real atmosphere, they will be smoothed out. Moreover, vacuum polarization can induce resonant mode conversion 
which may significantly affect the radiation transport (Lai & Ho 2001). We will present atmosphere models including vacuum 
polarization in a future work (Ho & Lai 2001). 

Finally, although in this paper we have included the medium effect on the polarization modes (see Section ^.f^ ), at 
sufficiently high densities or low photon energies, the refractive index of the medium deviates from unity or even becomes 
complex. This will introduce additional effects on the radiative transport and the resulting spectra. This issue will be addressed 
in a future work. 
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Lai 2001). 
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APPENDIX A: NORMAL MODE POLARIZATION VECTOR FOR MAGNETIZED ELECTRON-ION 
PLASMA 

We outline here the derivation of the polarization vector e-* for a magnetized electron-ion plasma. Many previous works on 
opacities did not include the ion effects on the polarization vector (e.g., Ventura 1979; Meszaros 1992; see, however, Bulik 
& Pavlov 1996). Our derivation is patterned after Shafranov (1967) (see also Ginzburg 1970), who calculated the index of 
refraction but did not give the explicit expression for the polarization vector. 

Neglecting the magnetic susceptibility and assuming no external charges or currents, the Maxwell equations are used to 
derive the equation for plane waves with E cx e'^'' ''""'-': 



[eij + n' {kikj — Ej 



0, 



(Al) 



where k = k/fc, Ei and ki are the ith component of the electric field and wave vector of the plane wave, respectively, and 
n = ck/ LJ is the refractive index. In the frame where the external magnetic field B is aligned along the z'-axis, the dielectric 
tensor is given by 



(A2) 



e 


ig 







e 











n 



where 



s 



^ Ai 



(A3) 
(A4) 

(A5) 



Here the sums run over each charged particle species s in the plasma, and u — (jjg^/uj^ and v = ajp^/tj^, where lobs and ujps 
are the cyclotron and plasma frequencies of charged particle s, respectively. Damping of the particle motion is accounted for 
in As = 1 + iVsjoj, where is the damping rate of a particle of species s. As in previous studies, damping is usually taken to 
be small, so that A^ 1. In the frame where the propagation wave vector k is aligned along the z-axis and B is in the x-z 
plane, the dielectric tensor is given by 



e cos^ 9b +1] sin^ 9b 
—ig cos 9b 



ig cos 9 b {e — rj) sm9B cos6b 
e —igsm9B 
{e — rj) sin 9 b cos 6 B igsinOs esiv? 9b + rj cos^ 9b 



(A6) 



Substituting equation (A6) into equation (Al) and setting the determinant of the resulting equation to zero yields the index 
of refraction. 



B±{B^ - AAC 



1/2 



2A 



where 

A = 
B = 
C = 



e sin^ 9b + Tj cos^ 6b 



g -ej7jsm 

2\ 



2eri 



g )v- 



(A7) 

(A8) 
(A9) 
(AlO) 



To obtain the polarization vector, we write the electric field as E = e-* = eo{iKj,l,iKz,j), where iKj = Ex/Ey 



iK, 



Ez/Ey^ and eo = (1 + K'^ + ^z,j) ^^"^ from the normalization |Ep = 1. From equation (Al), we find 



(All) 



which gives equation (2.44) when considering the case of electrons and one species of ions. Equation (Al) can also be used to 
obtain Kj, which is given by equation (2.43), with 



erj) sm 



2gri cos 9b 



(A12) 



Specializing to the case of electrons and one species of ions, equation ( |A12|) can be simplifi ed to equa tion ( 2.45 ). Once e-* is 
rotated into the frame with B along the z-axis, the polarization vector becomes equations (2.41) and (2.42). 
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Figure 1. Angle-averaged absorption opacities K^^^ (upper panels) and scattering opacities K^'^ (lower panels) as functions of energy for 
magnetic fields B = 4.7 X 10^^ G (left panels) and 10^^ G (right panels). The densities and temperatures are chosen to be representative 
of the hydrogen atmospheres (Z = A = 1) studied in Section W (see Fig. Bl): for B = 4.7 X 10^'^ G, Thomson depth r = 1 corresponds to 
p = 1 g/cm3 and T = 3.7 X lO" K and r = 22 corresponds to p = 13 g/cri? and T = 5.9 X 10^ K; for B = 10^^ G, Thomson depth t = 1 
corresponds to p = 1.1 g/cm^ and T = 3.4 X 10^ K and t = 680 corresponds to p = 370 g/cm^ and T = 7.0 X 10® K. The solid lines are 
for the extraordinary mode (j = 1) and the dotted lines are for the ordinary mode {j = 2). We note that the nonmagnetic opacities are 
similar to the ordinary mode opacities in t erms of energy dependence, and the variations in the ordinary mode opacities at low energies 
{E ^ ?iu)p) are due to K^j [oc lo^; see eq. ( ^.44| )]. 
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Figure 2. Convergence tests after ~20 global iterations (temperature corrections) for a fully ionized hydrogen atmosphere model with 
Tcft = 5 X 10® K, B = 10^ '^ G. and &b = 0. Fractional temperature correction (upper panel), deviation from radiative equilibrium 
[middle panel; see eqs. ( |3.3| )-( ^!^ )] , and deviation from constant total flux (lower panel) are plotted as functions of Thomson depth r. 
The larger deviations around r ~ 1 — 100 in the lower panel are due to the extraordinary-mode photons decoupling from the matter at 
T 100 (see Fig. a); this region has the slowest convergence, but the deviations can be reduced further with additional iterations. 
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Figure 3. Temperature T as a function of Tliomson deptli r (upper panel) and Rosseland mean depth t-r (lower panel) for fully ionized, 
pure hydrogen atmospheres with Toff = 5 X lO'' K. a nd mag netic fields B = 0, 4.7 X 10^^ iq14 q ^ q -pj^^ dashed lines show 

the grey temperature profile given by equation (B-U). 
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Figure 4. Temperature as a function of density for the various fully ionized atmosphere models considered in Section 1.3 The upper 
panel shows the H and He models at B = 0, 4.7 X 10^^ G, and 10^* G, all for T^g = 5 X lO'' K and 0^ = 0. The lower panel shows the 
H models at B = IQi"' G for T<,ff = 10*^ K, 5 X 10« K, and 7 X 10^ K. 
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Figure 5. The upper panels show the Thomson depth at which photons decouple from the matter [effective optical depth t^'J = 1; 
see eq. as a function of energy for the T^ff = 5 X 10'^ K, 0^ = 0, and B = 4.7 X lO^^ G (left panels) and lO^'* G (right panels) 

atmosphere models. Also plotted are the local temperature T (middle panels) and density p (lower panels) at these Thomson depths. 
The solid lines are for the extraordinary mode (j = 1), and the dashed lines are for the ordinary mode {j = 2). 
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Figure 6. Spectra of neutron star atmospheres with fully ionized hydrogen at T^g = 5 X 10^ K and ©s = 0. The dot-dashed line is 
for the B = 5 X 10^'^ G atmosphere, the solid line is for the B = 10^* G atmosphere, the short-dashed line is for the B = 4.7 X 10^^ G 
atmosphere, the long-dashed line is for the nonmagnetic atmosphere, and the dotted lines are for a blackbody with T = 5 X 10*^ K. At 
B = 5x 10" G, the Eni = 3.15 keV ion cyclotron feature (from ~2 to ~5 keV) has an equivalent width ~1.9 keV, while at _B = lO^^ G, 
the Egi = 0.63 keV feature (from ~0.4 to ~1 keV) has an equivalent width ~90 eV. Note that the extremely narrow natural/thermal 
lines at Esi have negligible equivalent widths. 
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Figure 7. Spectra of neutron star atmospheres with fully ionized hydrogen at T^g = 7 X 10*^ (upper panel) and lO'' K (lower panel) and 
©B = 0. The dot-dashed line is for the B = 5 X 10^* G atmosphere, the solid lines are for the B = 10^** G atmosphere, the short-dashed 
line is for the B = 4.7 X lO^'^ G atmosphere, the long-dashed lines are for the nonmagnetic atmosphere, and the dotted lines are for a 
blackbody with T = 7 X 10*^ K and 10*^ K. Note that the ion cyclotron resonance features are more pronounced when Esi ^ SfcTeff: at 
_B = 5 X IQl* G and Teff = 7 X 10'' K, the Esi = 3.15 keV feature (from ~2 to ~5 keV) has an equivalent width ~1.8 keV, while at 
B = 10^* G and T^g = 10^ K, the Esi = 0.63 keV feature (from ~0.4 to ~1 keV) has an equivalent width ~0.4 keV. Also note that the 
extremely narrow natural/thermal lines at Egi have negligible equivalent widths. 
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Figure 8. Spectra of neutron star atmospheres with fully ionized hydrogen at T^ff = 5 X 10 K and 0b = 0. The solid and dotted lines 
are for the B = 10^* G extraordinary mode (j = 1) and ordinary mode (j = 2), respectively, and the short and long-dashed lines are for 
the B = 4.7 X 10^'^ G extraordinary and ordinary modes, respectively. 
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Figure 9. Spectra of neutron star atmospheres with fully ionized helium at T^ff = 5 X 10^ K and @b = 0. The dot-dashed line is for 
the B = 5 X 10^'' G atmosphere, the solid line is for the B = 10^* G atmosphere, the short-dashed line is for the B = 4.7 X 10^^ G 
atmosphere, the long-dashed line is for the nonmagnetic atmosphere, and the dotted lines are for a blackbody with T = 5 X 10^ K. At 
B = 5x 10" G, the Esi = 1.58 keV ion cyclotron feature (from ~1 to ~2 keV) has an equivalent width ~0.51 keV, while at _B = 10^'' G, 
the Egi = 0.32 keV feature (from ~0.2 to ~0.4 keV) has an equivalent width '^10 cV. 
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Figure 10. Comparison of t he t emperature-density profiles between the full radiative transport and diffusion approximated transport 



models considered in Section 4.4 
the dotted line is for Qg = 90 



The solid line is for Qb 
(diffusion). 



: (full radiative transport), the dashed line is for = (diffusion), and 
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Figure 11. Spectra of neutron star atmospheres with fully ionized hydrogen at _B = 5 X 10^** G (upper panel) and 10^'* G (lower panel) 
and Tcft = 5 X 10*^ K. The solid lines are for Qb = (full radiative transport), the short-dashed lines are for @b = (diffusion), the 
dotted lines are for = 90° (diffusion), and the long-dashed lines are for a blackbody with 5 X 10® K. 



© 0000 RAS, MNRAS 000, 



